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Abstract. - We study topological aspects of a compact lattice superconductor, and show 
that the characteristic energy splitting, A, between almost degenerate ground states, is simply 
related to a novel order parameter W , which is closely related to large Wilson loops. Using 
Monte Carlo methods, we study the scaling properties of W close to the deconfining phase 
transition, and conclude that A ~ e~^^^, where L is the size of the system, thus giving 
quantitative support to the vortex tunneling scenario proposed by Wen. 



The concept of topological order was first introduced to describe the ground states of the 
quantum Hall system [1]. Although the present interest in topological order to large extent 
derives from the quest for exotic non- Fermi liquid states of relevance for the high- Tc supercon- 
ductors [2-6], the concept is of much more general interest. Both spin liquids [7,8], believed 
to be of relevance for frustrated magnets, and ordinary BCS superconductors with dynamical 
electromagnetism [9,10], are examples of topologically ordered systems. Two characteristics 
of topological order are particularly striking; fractionally charged quasiparticles, and a ground 
state degeneracy depending on the topology of the underlying space, which is lifted by tun- 
neling processes. In the case of the most celebrated example - the Laughlin FQH states - 
both these properties are well understood [11], but for superconductors, the situation is less 
clear. It was only in 1990 that Kivelson and Rokshar pointed out that the quasiparticles are 
indeed fractional [12], and to our knowledge there has been no quantitative study of ground 
state degeneracy and splitting. 

We improve on this situation by a numerical study of a 2D compact lattice superconductor, 
using a novel order parameter W , which is sensitive to vortex tunneling processes. Our results 
are in quantitative agreement with predictions based on topological order, and rules out the 
possibility of a spontaneously broken discrete symmetry. A compact lattice superconductor 
captures important properties of real type II superconductors [10], and the 2D case is also 
of intrinsic interest in the context of effective gauge theories for the high- Tc cuprates [2-6]. 
Examples of recent proposals are the U{1) slave-boson theory [2], the nodal liquid and the 
spinon-chargeon Z2 theory [3,4], and compact QED3 [5]. 
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We use the Villain form of the action [13], 

S=^Y1 (^A'^r + 2nnr^ - qArf^f + ^Yl (-^^^ + ^nkr^f , (1) 

where the charge q phase field 9^ is defined on the sites r — (xq, xi,X2), and the compact vector 
potential on the links (r,r + e^), of a Eucledian 3D cubic lattice. B^^ = ep,i,\S7 ^A^x is 
the dual field strength, which is defined on the links of the dual lattice (^). 

We shall in the following concentrate on g = 2 pertinent both to a Cooper paired BCS 
superconductor and to the more exotic pairing scenarios in effective gauge models [2,4-6], 
although we will keep a general q > 1 where appropriate. The phase diagram of the q — 2 
model is well established (^): In the strong coupling limit J — > oo, 1^ becomes a Z2 gauge 
theory, which has a phase transition at finite g. This transition between the "confined" phase 
at large g and the superconducting, or Higgs, phase at small g, can rigorously be shown to 
extend to finite J [14], and is believed to extend all the way to the 3D XY transition in the 
g limit. This is confirmed by recent simulations [15]. The line at J = corresponds to 
compact QED which is confined at all g because of instanton effects [16]. 

The topological structure of the theory, which is our present concern, is particularly simple 
in the g limit of the Z2 line, where the ground state degeneracy on, say, a torus can be 
understood in terms of "visons" describing Z2 magnetic fluxes through the "holes" [4,7]. At 
finite g the ground states are connected by processes where these vortices tunnel around the 
closed cycles of the torus, leading to a predicted energy splitting A ^ e~'^^l^ [9]. Close 
to the phase transition, the vortices begin to proliferate and this simple formula might well 
break down. Note that in the case of spontaneous breakdown of a discrete symmetry, there 
would also be tunneling processes connecting different ground states, but with a splitting 
A ~ e""^^ 1^" characteristic of domain wall tunneling [9]. 

To move away from the limiting cases, and to study the region close to the phase transition, 
which is of particular interest in the models of high- Tc superconductivity referred to above, 
numerical simulations are mandatory. It will be advantageous to use a dual formulation in 
terms of flux tubes and monopoles (i.e. instantons) [13,17] and express the partition function 
for |(TJ as Z = Z]{m,,Af,} '^^t^ 

5 27rV^ mr • Kr'nir' + g^A^iVrKr'A^r' ■ (2) 

Here nir € 1? is the vorticity on the dual links and A'r € Z is the number of monopoles 
on the dual sites. The flux lines are constrained to start or end on the monopoles only in 
quanta of g, i.e., V • nir = qA^r- The interaction is given by 14 = 14;e*'' '', with 

— X^iy 4 sin^(fc,y/2) + A^^, which for large distances is a screened coulomb interaction 
V{r) = e~''/^/47rr with a screening length A given by A~^ = Jgq^ . In the limit A — > the 
interaction reduces to an on-site interaction. 




where the flux quantum $0 = — • Since each +(— ) monopole has precisely q outgoing (incom- 
ing) flux lines, there are (for q> V) two distinct phases: In the superconducting phase, where 

(^)We use a continuum notation whenever convenient, and the lattice version should be obvious. 
(^)The phase diagram in the <? = 1 case is controversial. 
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Fig. 1 - Eucledian space-time configuration corresponding to a vortex tunneling around a closed cycle 
in the X2 direction on the torus. The operator W = (exp(i A ■ dr)) measures the flux through the 
loop surrounding the shaded area in the figure. The loop may be decomposed into two vertical 
Polyakov loops, Po and Pq, which cancel each other for periodic boundary conditions in the spatial 
directions, and two horizontal ones, Pi and P/, that differ only on a single twisted link belonging to 
the indicated column that is twisted by the gauge transformation ilio- 

flux lines cost a lot of energy, the monopoles are confined in neutral pairs bound together by 
q flux lines. In the other phase, the flux lines condense, connect many monopoles, and form 
a large connected tangle which percolates through the whole system. The electric properties 
of these phases are encoded in the behavior of a large Wilson loop W = (exp (i <f^ A ■ dr)), 
for fractionally charged test particles. In the percolating flux phase, where large flux loops 
dominate, the flux on any patch the size of a correlation length in a surface spanned by C, 
can be considered as a random variable, thus giving an area law, while in the superconducting 
phase, the small loops will only contribute close to the edge of C, giving a perimeter law [13]. 

Using periodic boundary conditions on the m:s, the total flux for a conflguration is just 
the sum of the m:s through a full cross-section S of the system in the direction /i, ^oMf^ = 
$0 X]re5 "^r/i- Clearly ~ ^g**"*^'') is closely related to large Wilson loops and directly 
measures the presence of percolating flux, so we expect it to have area law behaviour in 
the confined phase. In the superconducting phase, however, W and W differ - the latter is 
insensitive to small flux loops and should not obey a perimeter law. For these reasons, we 
submit that is a good (non-local) order parameter for the deconfinement transition, and if 
the transition is continuous, we expect a finite size scaling relation, 

iy^(5,L) = (e'^^''*") = IT±(L/0, ^~|<Sr^ (4) 

to hold. Here 5 is a tuning parameter of the transition, e.g., S = (J — Jc)/Jc, W± is a 
universal scaling function, and ^ is the correlation length which diverges as the transition is 
approached (■^). 

Physically, configurations having non-zero total flux $o^/i around spatial directions of 
the torus correspond to vortex tunneling events. In order to include these in the partition 



(^)For q = 1 this would not give any information, since in this case W = 1 identically. 
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function we must be careful about boundary conditions. Returning to the original formulation 
in terms of A and 9, we note that periodic boundary conditions in the /i and v directions are 
compatible only with the total vorticity in the iiv plane being an integer multiple of q. We 
can, however, choose twisted boundary conditions, 

where fi^j/ is a large gauge transformation, that shifts on all links in one particular column 
in the a direction, {^va cyclic) by 2Tmf^i,/q, with n^j, = 1,. . . ,q. The twisted links corre- 
sponding to 17 10 are shown in Fig.^ Acting as an operator on states in the Hilbert space ftfj.i, 
changes the value of the Polyakov loop operator ~ exp(i ^ A^) (the integration is around a 
closed cycle on the torus), according to the topological algebra ilyf^Pi, — e'^'^"'"'/'?P,yri,y^ — 0. 
With the boundary conditions ijSJ), the total vorticity in the a direction is n^i, + n^f^ (mod q). 

In order to discuss quantum mechanical states, we now distinguish between the Euclidean 
time direction /i = 0, and the spatial directions i = 1,2. For simplicity we also specialize to 
q — 2 and put Qtj ~ Periodicity in space implies even total flux in the 0-direction, 

and, absent vortices, the Polyakov loop operators. Pi can be diagonalized, with eigenvalues 
±1, so the ground state manifold is spanned by the states |/i, /2), where ftio |1, 1) = |0, 1) etc. 
Tunneling matrix elements between these states can be extracted from the twisted partition 
functions, 

ZIm = Tr [e-^"n,o] . (6) 

The corresponding Euclidean path integral with the boundary condition (jSJ in the time direc- 
tion will, in the flux-monopole formulation, amount to summing over configurations with odd 
total flux through the (Oi) plane. Similarly, the untwisted partition function, Zl^ , corresponds 
to even total flux. Fig. ^ illustrates that Z^^^^ indeed corresponds to vortex-antivortex pairs 
tunneling in the 2-direction. In the presence of vortex loops and monopole-antimonopole pairs, 
the ground states are no longer eigenstates of the loop operators Pi{xo, xj). Nevertheless, by 
continuity, the degenerate ground states will persist, and since the effects of flux loops and 
monopoles in the superconducting phase are local, the tunneling between them is still given 
solely by the flux lines traversing the whole torus. In a real superconductor with unit charge 
quasiparticles present, the situation is more complicated, and the ground state is determined 
by competing tunneling processes [10]. 

Allowing, for simplicity, tunneling only in the 2-direction and restricting the Hilbert space 
to the two ground states (which is appropriate at low temperatures), we may write down a 
tunneling Hamiltonian with elements Hn = H22 — Eq, H12 = H21 = —A in the flux basis 
1/1), which is thus diagonal in the eigenstates |±) — (|0) ± |l))/\/2 of flw, with eigenvalues 
E± = Ea ^ A. From the deflnition of the twisted and untwisted partition function, it now 
follows, 

(±1 e-^" |±) = ± Zodd = e-''(^''T^) , (7) 
which gives the flnal formula for the energy splitting. 



_2/3A Zf^y Zofifi 



Zp,,_ + Zr 



dd 



= , (8) 



where the last equality follows since the total flux operator ^qM^ takes the values and tt 
for q = 2. Not only does this result allow the gap A to be explicitly calculated from W . 
It also implies, via Eq. Q), a scaling relation, A = A±{L/£^)/L. Expressed in the original 
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Fig. 2 - W as function of g for the A — > case (a), and of J'^ for A = 0.5 (b). The different curves 
cross right at the phase transition, where W becomes independent of the system size, according to 
Eq. @. 



variables, W takes the form W = (^Pi{0,X2)P{iP,X2)j, c.f. Fig.[T](4). The generalization to 
the full 4 state system is straightforward - the splitting pattern is now (— 2A, 0, 0, 2A), while 
the relation ||SJ) between A and W still holds. It is also possible to generalize to higher q > 1 
and other topologies. 

Deep into the superconducting phase we can estimate W by assuming that straight, and 
statistically independent, paths dominate the partition sum. An easy calculation using ^ 
gives Aq = [L/ a?)e^'^^^°- , where a is the lattice spacing, and c is the line tension of the flux 
lines (c — <&o/2.9 for A ^ a and j7rln(A/a) for A ^ a). To extend this analysis away from 
the small A and small g part of the phase diagram, and especially to the region close to the 
phase transition, we match the splitting Aq deep in the phase to the scaling expression given 
in Eq. Q), leading to 

A = |e-'V« , (9) 

for ^ < i, where ^ ^ \^\~'^ is now a physical correlation length. Thus the energy splitting gets 
renormalized and eventually closes as the transition is approached. The exponentially small 
splitting has been predicted by Wen, whereas here we obtain also a nontrivial prefactor. 

Equation ^ enables us to explicitly calculate the splitting A numerically using Monte 
Carlo simulations. We simulate the system in the flux line - monopole representation given 
by Eq. ^ or (|3Jl, using a variant of a recently described worm cluster update Monte Carlo 
algorithm [18], properly adapted to long range interactions and the existence of magnetic 
monopoles. This algorithm naturally includes global moves which change (i.e., we allow 
twists in the boundary conditions that are necessary for the calculation of W, whereas in a 
conventional Metropolis algorithm the acceptance ratio for such moves becomes exponentially 
small with increasing lattice size. For convenience we allow twists not only in the time direction 
but also in space, i.e., we simulate the system in a grand canonical ensemble. We have checked 
that this does not change our results. The details of the numerical methods will be described 
elsewhere [19]. 

We have carried out simulations for constant A 0, 0.5, 1, 2, varying J in each simulation. 
For a given A, the critical value Jc (or in case A = 0) of the deconfinemcnt transition is 



('')This operator differs from 1 only at the one particular link where the A field is twisted. This explains, in 
terms of the original variables, why W has no perimeter law. 
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Fig. 3 - Test of Eq. Q for A = 0. When plotted against L/^, the combination 2A^^/L falls onto 
two separate branches. The data for g < Qc shows an exponential dependence (indicated by the solid 
straight line) in agreement with Eq. 0. The data for g > gc shows an area law indicated by the 
horizontal dotted line. Black symbols are for (3 — L and red (gray) symbols are for /3 = 2L. 



located using the nonlocal order parameter W, as shown in Fig. |21 According to Eq. W 
for different system sizes should cross right at the transition. In the topologically ordered 
phase W tends to one, implying that A — > and the states becomes close to degenerate. The 
exponential dependence of A on L in Eq. ((HJ is obtained in Fig. |31 for A = and in Fig. 01 
for A = 0.5. We plot 2A^^/L vs L/£^, where the correlation length is given by ^ = 
The correlation length exponent i/ is adjusted until the data collapse onto two branches (for 
positive and negative S), which happens for v = 0.63, consistent with the 3D Ising critical 
behavior expected for g = 2. The proportionality constant A is fixed by Eq. 1^. The straight 




Fig. 4 - Same as Fig. 01 but for A = 0.5. 
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line in the lin-log plot indeed demonstrates the exponential scaling in the topologically ordered 
phase, whereas an area law indicated by the horizontal line is obtained in the confined phase. 
Note that a e~'^'^ behavior characteristic of a spontaneously broken discrete symmetry 
is excluded. The results for A = 1 and 2 are very similar. In principle, the simulations 
done for {3 = L cannot rule out a weak dependence on j3/L in Eq. (jHl- Therefore we have 
also calculated A from anisotropic systems with /? = 2L, see Fig. showing that any such 
residual dependence on (3/ L must be extremely weak. 

In summary, we have shown that the ground state energy splitting A, between the almost 
degenerate ground states in the topologically ordered phase of a lattice superconductor is 
simply related to a novel nonlocal order parameter W of the deconfining transition. Using a 
Monte Carlo algorithm we calculated W , demonstrated that it is a good order parameter for 
the phase transition, and established that A has an exponential dependence on system size in 
full agreement with the vortex tunneling mechanism proposed by Wen. This suggests that the 
lattice superconductor shows the characteristics of topological order all the way to the phase 
transition. Let us finally remark that the methods used in this paper should be applicable 
also to a model with quenched impurities. 

* * * 
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